// Copyright (c) Lawrence Livermore National Security, LLC and other VisIt
// Project developers.  See the top-level LICENSE file for dates and other
// details.  No copyright assignment is required to contribute to VisIt.

// ************************************************************************* //
//                            avtvolimageFileFormat.C                        //
// ************************************************************************* //

#include <avtvolimageFileFormat.h>

#include <string>

#include <vtkFloatArray.h>
#include <vtkRectilinearGrid.h>
#include <vtkStructuredGrid.h>
#include <vtkUnstructuredGrid.h>

#include <avtDatabaseMetaData.h>

#include <Expression.h>

#include <InvalidVariableException.h>
#include <InvalidDBTypeException.h>

#include <DebugStream.h>

#include <fcntl.h>
#include <visit-config.h>

using namespace std;

#ifdef _WIN32
#include <io.h>
#define OPEN  ::_open
#define CLOSE ::_close
#define LSEEK ::_lseek
#else
#include <unistd.h>
#define O_BINARY 0
#define OPEN  ::open
#define CLOSE ::close
#define LSEEK ::lseek
#endif


// ****************************************************************************
//  Method: avtvolimageFileFormat constructor
//
//  Programmer: sjogreen2 -- generated by xml2avt
//  Creation:   Mon Jun 20 08:59:36 PDT 2011
//
// ****************************************************************************

avtvolimageFileFormat::avtvolimageFileFormat(const char *filename)
    : avtSTMDFileFormat(&filename, 1)
{
    // INITIALIZE DATA MEMBERS
    m_filename = filename;

// gridfilename is read from the header in the solution file

    m_initialized = false;
    m_data_stored = false; // set to true if you want the solution field to be stored internally in this object
    m_grid_stored = false; // not currently used
}


// ****************************************************************************
//  Method: avtvolimageFileFormat::FreeUpResources
//
//  Purpose:
//      When VisIt is done focusing on a particular timestep, it asks that
//      timestep to free up any resources (memory, file descriptors) that
//      it has associated with it.  This method is the mechanism for doing
//      that.
//
//  Programmer: sjogreen2 -- generated by xml2avt
//  Creation:   Mon Jun 20 08:59:36 PDT 2011
//
// ****************************************************************************

void
avtvolimageFileFormat::FreeUpResources(void)
{
}


// ****************************************************************************
//  Method: avtvolimageFileFormat::PopulateDatabaseMetaData
//
//  Purpose:
//      This database meta-data object is like a table of contents for the
//      file.  By populating it, you are telling the rest of VisIt what
//      information it can request from you.
//
//  Programmer: sjogreen2 -- generated by xml2avt
//  Creation:   Mon Jun 20 08:59:36 PDT 2011
//
// ****************************************************************************

void
avtvolimageFileFormat::PopulateDatabaseMetaData(avtDatabaseMetaData *md)
{
    debug5 << "Populatedatabase begins" << endl;

    if( !m_initialized )
        Initialize();

    // Mesh information
    avtMeshMetaData *mmd = new avtMeshMetaData;

    if (m_CartGrid)
    {
        mmd->name = "Cartesian grid";
    }
    else
    {
        mmd->name = m_gridname;
    }

    mmd->spatialDimension = 3;
    mmd->topologicalDimension = 3;
    mmd->numBlocks = m_nblocks;
    //   mmd->hasUnits = true;
    //   mmd->units = "meters";

    if (m_CartGrid)
    {
        mmd->meshType = AVT_RECTILINEAR_MESH;
    }
    else
    {
        mmd->meshType = AVT_CURVILINEAR_MESH;
    }

    debug5 << "before b-loop " << endl;
    mmd->blockNames.resize(m_nblocks);
    for( int b=0 ; b < m_nblocks ; b++ )
    {
        debug5 << "Block = " << b << endl;
        debug5 << "grid size = " << m_gridsize[b] << endl;
        debug5 << "no of pts = " << m_ni[b] << " " << m_nj[b] << " " << m_nk[b] << endl;
        char buf[50];
        snprintf(buf,50,"h=%f",m_gridsize[b]);
        string bname= buf;
        mmd->blockNames[b] = bname;
    }
    md->Add(mmd);
    debug5 << "after b-loop " << endl;

    // Scalar variable information
    avtScalarMetaData *smd1 = new avtScalarMetaData;
    smd1->name = m_mode;
// copy the mesh name from above
    smd1->meshName = mmd->name;
    smd1->centering = AVT_NODECENT;
    smd1->hasUnits = false;
    md->Add(smd1);
    debug5 << "Populatedatabase done " << endl;
}


// ****************************************************************************
//  Method: avtvolimageFileFormat::GetMesh
//
//  Purpose:
//      Gets the mesh associated with this file.  The mesh is returned as a
//      derived type of vtkDataSet (ie vtkRectilinearGrid, vtkStructuredGrid,
//      vtkUnstructuredGrid, etc).
//
//  Arguments:
//      domain      The index of the domain.  If there are NDomains, this
//                  value is guaranteed to be between 0 and NDomains-1,
//                  regardless of block origin.
//      meshname    The name of the mesh of interest.  This can be ignored if
//                  there is only one mesh.
//
//  Programmer: sjogreen2 -- generated by xml2avt
//  Creation:   Mon Jun 20 08:59:36 PDT 2011
//
//  Modifications:
//    Brad Whitlock, Mon Aug 15 09:33:38 PDT 2011
//    Free resources before throwing exceptions.
//
// ****************************************************************************

vtkDataSet *
avtvolimageFileFormat::GetMesh(int domain, const char *meshname)
{
//
// **************NEED TO CHANGE TO A STRUCTURED MESH**************
//
//  COULD ENCODE THE NAME OF THE GRID FILE INTO THE HEADER OF THE SOLUTION FILE AND
//  PASS THIS INFO THROUGH MESHNAME
//
    if( !m_initialized )
        Initialize();
    //    YOU MUST IMPLEMENT THIS
    debug5 << "entering get mesh " << endl;
    debug5 << "meshname = '" << meshname << "'" << " domain = " << domain << endl;

    if (m_CartGrid)
    {
        debug5 << "constructing a Cartesian grid" << endl;

        vtkFloatArray *coords[3]= {0,0,0};
        coords[0] = vtkFloatArray::New();
        coords[0]->SetNumberOfTuples(m_ni[domain]);
        coords[1] = vtkFloatArray::New();
        coords[1]->SetNumberOfTuples(m_nj[domain]);
        coords[2] = vtkFloatArray::New();
        coords[2]->SetNumberOfTuples(m_nk[domain]);
        float* x = (float *)coords[0]->GetVoidPointer(0);
        float* y = (float *)coords[1]->GetVoidPointer(0);
        float* z = (float *)coords[2]->GetVoidPointer(0);
        for( int i=0 ; i < m_ni[domain] ; i++ )
            x[i] = i*m_gridsize[domain]+m_xmin[domain];
        for( int i=0 ; i < m_nj[domain] ; i++ )
            y[i] = i*m_gridsize[domain]+m_ymin[domain];
        for( int i=0 ; i < m_nk[domain] ; i++ )
            z[i] = i*m_gridsize[domain]+m_zmin[domain];
        vtkRectilinearGrid *rgrid= vtkRectilinearGrid::New();
        int dims[3];
        dims[0] = m_ni[domain];
        dims[1] = m_nj[domain];
        dims[2] = m_nk[domain];
        rgrid->SetDimensions(dims);
        rgrid->SetXCoordinates(coords[0]);
        rgrid->SetYCoordinates(coords[1]);
        rgrid->SetZCoordinates(coords[2]);
        coords[0]->Delete();
        coords[1]->Delete();
        coords[2]->Delete();
        debug5 << "done get mesh " << endl;
        return rgrid;
    }
    else // curvilinear grid in Z
    {
        debug5 << "constructing a curvilinear grid" << endl;

        //int ndims=3;
        int dims[3] = {1,1,1};
        dims[0] = m_ni[domain];
        dims[1] = m_nj[domain];
        dims[2] = m_nk[domain];
        int nnodes = dims[0]*dims[1]*dims[2];
        debug5 << "dims[0] = " << dims[0] << " dims[1] = " << dims[1] << " dims[2] = " << dims[2] << endl;

        float *zarray = new float[nnodes];

// POPULATE THE zarray variable with grid point coordinates...
        //int pnt=0;

//      for(int k = 0; k < dims[2]; ++k)
//        for(int j = 0; j < dims[1]; ++j)
//      for(int i = 0; i < dims[0]; ++i)
//      {
// // cartesian mesh (for testing)
//        zarray[pnt] =  k*m_gridsize[domain]+m_zmin[domain];
//        pnt++;
//      }

// read z-coordinate from grid file
        char errmsg[500];
        int fd = OPEN(m_gridfilename.c_str(),O_RDONLY|O_BINARY);
        if( fd == -1 )
        {
            delete [] zarray;
            snprintf(errmsg,500,"Error opening grid file %s",m_gridfilename.c_str());
            EXCEPTION1( InvalidDBTypeException, errmsg );
        }

// get the data
        off_t nr = LSEEK(fd,m_gridoffset[domain],SEEK_CUR);
        if( nr != m_gridoffset[domain] )
        {
            CLOSE(fd);
            delete [] zarray;
            snprintf(errmsg,500,"Error accessing array in %s",m_gridfilename.c_str());
            EXCEPTION1( InvalidDBTypeException, errmsg );
        }
        if( m_prec == 4 )
        {
            nr = read(fd,zarray,sizeof(float)*nnodes); // read straight into zarray
            if( (size_t)nr != sizeof(float)*nnodes )
            {
                CLOSE(fd);
                delete [] zarray;
                snprintf(errmsg,500,"Error reading grid array in %s" , m_gridfilename.c_str());
                EXCEPTION1( InvalidDBTypeException, errmsg );
            }
        }
        else
        {
            double* tmp=new double[nnodes];
            nr = read(fd,tmp,sizeof(double)*nnodes); // read into tmp (double) array
            if( (size_t)nr != sizeof(double)*nnodes )
            {
                CLOSE(fd);
                delete [] zarray;
                delete [] tmp;
                snprintf(errmsg,500,"Error reading dp grid array in %s" , m_gridfilename.c_str());
                EXCEPTION1( InvalidDBTypeException, errmsg );
            }
            for( size_t i = 0 ; i < (size_t)nnodes ; i++ ) // copy over to zarray
                zarray[i] = tmp[i];
            delete[] tmp;
        }
// end reading z-ccordinates: close the file
        CLOSE(fd);

        vtkStructuredGrid *sgrid = vtkStructuredGrid::New();
        vtkPoints * points = vtkPoints::New();
        sgrid->SetPoints(points);
        sgrid->SetDimensions(dims);
        points->Delete();
        points->SetNumberOfPoints(nnodes);

        float *pts = (float *) points->GetVoidPointer(0);
        float *zc = zarray;

        for(int k = 0; k < dims[2]; ++k)
            for(int j = 0; j < dims[1]; ++j)
                for(int i = 0; i < dims[0]; ++i)
                {
                    *pts++ = i*m_gridsize[domain]+m_xmin[domain]; // Cartesian in x and y
                    *pts++ = j*m_gridsize[domain]+m_ymin[domain];
                    *pts++ = *zc++;
                }

// delete temporary arrays
        delete [] zarray;

        return sgrid;
    }

}


// ****************************************************************************
//  Method: avtvolimageFileFormat::GetVar
//
//  Purpose:
//      Gets a scalar variable associated with this file.  Although VTK has
//      support for many different types, the best bet is vtkFloatArray, since
//      that is supported everywhere through VisIt.
//
//  Arguments:
//      domain     The index of the domain.  If there are NDomains, this
//                 value is guaranteed to be between 0 and NDomains-1,
//                 regardless of block origin.
//      varname    The name of the variable requested.
//
//  Programmer: sjogreen2 -- generated by xml2avt
//  Creation:   Mon Jun 20 08:59:36 PDT 2011
//
//  Modifications:
//    Brad Whitlock, Mon Aug 15 09:31:54 PDT 2011
//    Free resources before throwing exceptions.
//
// ****************************************************************************

vtkDataArray *
avtvolimageFileFormat::GetVar(int domain, const char *varname)
{
    if( !m_initialized )
        Initialize();
    debug5 << "entering get var " << endl;
    char errmsg[500];

    //    YOU MUST IMPLEMENT THIS
    vtkFloatArray* arr = vtkFloatArray::New();
    size_t npts = ((size_t) m_ni[domain])*m_nj[domain]*m_nk[domain];
    arr->SetNumberOfTuples(npts);
    float* data = (float *)arr->GetVoidPointer(0);
    if( m_data_stored )
    {
        float *ptr = m_dataptr[domain];
        for( size_t i = 0 ; i < npts ; i++ )
            data[i] = ptr[i];
    }
    else
    {
        int fd = OPEN(m_filename.c_str(),O_RDONLY|O_BINARY);
        if( fd == -1 )
        {
            snprintf(errmsg,500,"Error opening file %s",m_filename.c_str());
            EXCEPTION1( InvalidDBTypeException, errmsg );
        }
        off_t nr = LSEEK(fd,m_offset[domain],SEEK_CUR);
        if( nr != m_offset[domain] )
        {
            CLOSE(fd);
            snprintf(errmsg,500,"Error accessing array in %s",m_filename.c_str());
            EXCEPTION1( InvalidDBTypeException, errmsg );
        }
        if( m_prec == 4 )
        {
            nr = read(fd,data,sizeof(float)*npts);
            if( (size_t)nr != sizeof(float)*npts )
            {
                CLOSE(fd);
                snprintf(errmsg,500,"Error reading array in %s" , m_filename.c_str());
                EXCEPTION1( InvalidDBTypeException, errmsg );
            }
        }
        else
        {
            double* tmp=new double[npts];
            nr = read(fd,tmp,sizeof(double)*npts);
            if( (size_t)nr != sizeof(double)*npts )
            {
                CLOSE(fd); 
                delete [] tmp;
                snprintf(errmsg,500,"Error reading array in %s" , m_filename.c_str());
                EXCEPTION1( InvalidDBTypeException, errmsg );
            }
            for( size_t i = 0 ; i < npts ; i++ )
                data[i] = tmp[i];
            delete[] tmp;
        }
// close the file
        CLOSE(fd);

    }
    debug5 << "done get var " << endl;
    return arr;
}


// ****************************************************************************
//  Method: avtvolimageFileFormat::GetVectorVar
//
//  Purpose:
//      Gets a vector variable associated with this file.  Although VTK has
//      support for many different types, the best bet is vtkFloatArray, since
//      that is supported everywhere through VisIt.
//
//  Arguments:
//      domain     The index of the domain.  If there are NDomains, this
//                 value is guaranteed to be between 0 and NDomains-1,
//                 regardless of block origin.
//      varname    The name of the variable requested.
//
//  Programmer: sjogreen2 -- generated by xml2avt
//  Creation:   Mon Jun 20 08:59:36 PDT 2011
//
//  Modifications:
//
// ****************************************************************************

vtkDataArray *
avtvolimageFileFormat::GetVectorVar(int domain, const char *varname)
{
    return 0;
}

// ****************************************************************************
//  Method: avtvolimageFileFormat::Initialize
//
//  Purpose:
//      Read the file.
//
//  Programmer: sjogreen2 -- generated by xml2avt
//  Creation:   Mon Jun 20 08:59:36 PDT 2011
//
//  Modifications:
//    Brad Whitlock, Fri Aug 12 11:49:32 PDT 2011
//    Prepend the path to the file to the grid file so we can find it.
//
// ****************************************************************************

void avtvolimageFileFormat::Initialize()
{
    char errmsg[500];
    debug5 << "Initialize begins" << endl;
    size_t k=m_filename.rfind('.');
    if( k != string::npos )
        m_mode = m_filename.substr(k+1);
    else
        m_mode = " ";
    if( !(m_mode == "ux" || m_mode == "uy" || m_mode == "uz" || m_mode == "rho" ||
            m_mode == "p" || m_mode == "s" || m_mode == "div" || m_mode == "curl" ||
            m_mode == "mag" || m_mode == "veldiv" || m_mode == "velcurl" ||
            m_mode == "velmag") ) // AP removed .z extension
    {
        snprintf(errmsg,500,"Error: Unknown volimage mode %s" , m_mode.c_str() );
        EXCEPTION1( InvalidDBTypeException, errmsg );
    }
    debug5 << "mode = " << m_mode << endl;

    int fd = OPEN( m_filename.c_str(), O_RDONLY|O_BINARY );
    if( fd == -1 )
    {
        snprintf(errmsg,500,"Error opening file %s",m_filename.c_str());
        EXCEPTION1( InvalidDBTypeException, errmsg );
    }
    debug5 << "file opened " << endl;
    size_t nr = read(fd,&m_prec,sizeof(int) );
    if( nr != sizeof(int) )
    {
        snprintf(errmsg,500,"Error reading precision in %s",m_filename.c_str());
        EXCEPTION1( InvalidDBTypeException, errmsg );
    }
    nr = read(fd,&m_nblocks,sizeof(int) );
    if( nr != sizeof(int) )
    {
        snprintf(errmsg,500,"Error reading nblocks in %s",m_filename.c_str());
        EXCEPTION1( InvalidDBTypeException, errmsg );
    }
    nr = read(fd,&m_time,sizeof(double) );
    if( nr != sizeof(double) )
    {
        snprintf(errmsg,500,"Error reading time in %s",m_filename.c_str());
        EXCEPTION1( InvalidDBTypeException, errmsg );
    }
    debug5 << "prec = " << m_prec << " " << " nblocks = " << m_nblocks << " time= " << m_time << endl;

// read number of characters in gridfile name
    int nchars;
    char cgfn[200];
    nr = read(fd,&nchars,sizeof(int) );
    if( nr != sizeof(int) )
    {
        snprintf(errmsg,500,"Error reading nchars in %s",m_filename.c_str());
        EXCEPTION1( InvalidDBTypeException, errmsg );
    }

// was there a grid file?
    m_CartGrid = (nchars == 0);

    if (!m_CartGrid)
    {
// read name of gridfile
        nr = read(fd,&cgfn,nchars*sizeof(char) );
        if( nr != nchars*sizeof(char) )
        {
            snprintf(errmsg,500,"Error reading grid file name in %s",m_filename.c_str());
            EXCEPTION1( InvalidDBTypeException, errmsg );
        }
        cgfn[nchars]='\0';
        debug5 << "reading grid file name nchars= " << nchars << " string='" << cgfn << "'" << endl;

        // Get the path to the filename that we're opening.
        std::string path(m_filename);
        std::string::size_type pos = path.rfind(VISIT_SLASH_STRING);
        if(pos != std::string::npos)
        {
            path = path.substr(0, pos);
            if(path.size() > 0 && path[path.size()-1] != VISIT_SLASH_CHAR)
                path += VISIT_SLASH_STRING;
        }
        else
            path = "";

#ifndef _WIN32
        // If the first character of the string is slash then don't prepend path.
        if(cgfn[0] == VISIT_SLASH_CHAR)
            path = "";
#endif

// grid file name
        m_gridname = cgfn;
        m_gridfilename = path + cgfn;
        debug5 << "m_gridfilename = '" << m_gridfilename << "'" << endl;
    }
    else
    {
        debug5 << "No grid file, assuming Cartesian grid" << endl;
    }

    m_gridsize.resize(m_nblocks);
    m_xmin.resize(m_nblocks);
    m_ymin.resize(m_nblocks);
    m_zmin.resize(m_nblocks);
    m_ni.resize(m_nblocks);
    m_nj.resize(m_nblocks);
    m_nk.resize(m_nblocks);
    m_offset.resize(m_nblocks);
    m_gridoffset.resize(m_nblocks);
    off_t header_offset = 3*sizeof(int)+sizeof(double)+nchars*sizeof(char);
    int dims[6];
    debug5 << "reading block headers " << endl;
    for( int b=0 ; b < m_nblocks ; b++ )
    {
        debug5 << "b = " << b << endl;
        nr = read(fd,&m_gridsize[b],sizeof(double));
        if( nr != sizeof(double) )
        {
            snprintf(errmsg,500,"Error reading gridsizes in %s",m_filename.c_str());
            EXCEPTION1( InvalidDBTypeException, errmsg );
        }
        nr = read(fd,&m_xmin[b],sizeof(double));
        if( nr != sizeof(double) )
        {
            snprintf(errmsg,500,"Error reading xmin in %s",m_filename.c_str());
            EXCEPTION1( InvalidDBTypeException, errmsg );
        }
        nr = read(fd,&m_ymin[b],sizeof(double));
        if( nr != sizeof(double) )
        {
            snprintf(errmsg,500,"Error reading ymin in %s",m_filename.c_str());
            EXCEPTION1( InvalidDBTypeException, errmsg );
        }
        nr = read(fd,&m_zmin[b],sizeof(double));
        if( nr != sizeof(double) )
        {
            snprintf(errmsg,500,"Error reading zmin in %s",m_filename.c_str());
            EXCEPTION1( InvalidDBTypeException, errmsg );
        }
        nr = read(fd,dims,sizeof(int)*6);
        if( nr != sizeof(int)*6 )
        {
            snprintf(errmsg,500,"Error reading dimensions in %s",m_filename.c_str());
            EXCEPTION1( InvalidDBTypeException, errmsg );
        }
        m_ni[b] = dims[1]-dims[0]+1;
        m_nj[b] = dims[3]-dims[2]+1;
        m_nk[b] = dims[5]-dims[4]+1;
        m_xmin[b] += (dims[0]-1)*m_gridsize[b];
        m_ymin[b] += (dims[2]-1)*m_gridsize[b];
        m_zmin[b] += (dims[4]-1)*m_gridsize[b];
        header_offset += 4*sizeof(double) + 6*sizeof(int);
        debug5 << "b = " << b << " dims[0] = " << dims[0] << " dims[1] = " << dims[1] << " dims[2] = " << dims[2] <<
               " dims[3] = " << dims[3] << " dims[4] = " << dims[4] << " dims[5] = " << dims[5] << endl;
        debug5 << "x_min = " << m_xmin[b] << " y_min = " << m_ymin[b] << " z_min = " << m_zmin[b] << endl;
    }
    m_offset[0] = header_offset;
    int datasize;
    if( m_prec == 4 )
        datasize = sizeof(float);
    else
        datasize = sizeof(double);

// calculate offsets for accessing data
    for( int b=1 ; b < m_nblocks ; b++ )
        m_offset[b] = m_offset[b-1] + ((size_t) m_ni[b-1])*m_nj[b-1]*m_nk[b-1]*datasize;
    debug5 << " done headers, data stored =  " << m_data_stored << endl;

    if( m_data_stored )
    {
        m_dataptr.resize(m_nblocks);
        for( int b=0 ; b < m_nblocks ; b++ )
        {
            size_t npts = ((size_t) m_ni[b])*m_nj[b]*m_nk[b];
            m_dataptr[b] = new float[npts];
            if( m_prec == 4 )
            {
                nr = read(fd,&m_dataptr[b],npts*sizeof(float));
                if( nr != npts*sizeof(float) )
                {
                    snprintf(errmsg,500,"Error reading float data array %s",m_filename.c_str());
                    EXCEPTION1( InvalidDBTypeException, errmsg );
                }
            }
            else
            {
                double* tmp = new double[npts];
                nr = read(fd,&tmp,npts*sizeof(double));
                if( nr != npts*sizeof(double) )
                {
                    snprintf(errmsg,500,"Error reading double data array %s",m_filename.c_str());
                    EXCEPTION1( InvalidDBTypeException, errmsg );
                }
                for( size_t i = 0 ; i < npts ; i++ )
                    m_dataptr[b][i] = tmp[i];
                delete[] tmp;
            }
        }

    } // end data_stored

    CLOSE(fd); // closing the solution file

// now we need to open the gridfile and calculate its offsets,
// because the number of characters in the filename string might be different

    if (!m_CartGrid)
    {
        fd = OPEN(m_gridfilename.c_str(),O_RDONLY|O_BINARY);
        if( fd == -1 )
        {
            snprintf(errmsg,500,"Error opening grid file %s",m_gridfilename.c_str());
//       EXCEPTION1( InvalidDBTypeException, errmsg );
            m_CartGrid = true;
            debug5 << errmsg << endl;
            debug5 << "reverting to Cartesian Grid " << endl;
        }
    }

    if (!m_CartGrid)
    {
        debug5 << "gridfile opened " << endl;
        int dum_int;
        double dum_double;

        nr = read(fd,&dum_int,sizeof(int) );
        if( nr != sizeof(int) )
        {
            snprintf(errmsg,500,"Error reading precision in %s",m_gridfilename.c_str());
            EXCEPTION1( InvalidDBTypeException, errmsg );
        }
        nr = read(fd,&dum_int,sizeof(int) );
        if( nr != sizeof(int) )
        {
            snprintf(errmsg,500,"Error reading nblocks in %s",m_gridfilename.c_str());
            EXCEPTION1( InvalidDBTypeException, errmsg );
        }
        nr = read(fd,&dum_double,sizeof(double) );
        if( nr != sizeof(double) )
        {
            snprintf(errmsg,500,"Error reading time in %s",m_gridfilename.c_str());
            EXCEPTION1( InvalidDBTypeException, errmsg );
        }

// read number of characters in gridfile name
        nr = read(fd,&nchars,sizeof(int) );
        if( nr != sizeof(int) )
        {
            snprintf(errmsg,500,"Error reading nchars in %s",m_gridfilename.c_str());
            EXCEPTION1( InvalidDBTypeException, errmsg );
        }
        debug5 << "grid file nchars = " << nchars << endl;

        int grid_header_offset = 3*sizeof(int) + sizeof(double) + nchars*sizeof(char) +
                                 m_nblocks*(4*sizeof(double)+6*sizeof(int));

// calculate offsets for accessing data
        m_gridoffset[0] = grid_header_offset;

// we assume the same number of grid points in the grid file, and the same precision (float/double)
        for( int b=1 ; b < m_nblocks ; b++ )
            m_gridoffset[b] = m_gridoffset[b-1] + ((size_t) m_ni[b-1])*m_nj[b-1]*m_nk[b-1]*datasize;

        CLOSE(fd); // closing the grid file
    } // end if !m_CartGrid

//
    m_initialized = true;
    debug5 << "Initialize done " << endl;
}
